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ABSTRACT 

A three-step hybrid analysis technique, which successively uses the regular perturbation 
expansion method, the Pade expansion method, and then a Galerkin approximation, is 
presented and applied to some model boundary value problems. In the first step of the 
method, the regular perturbation method is used to construct an approximation to the 
solution in the form of a finite power series in a small parameter e associated with the 
problem. In the second step of the method, the series approximation obtained in step one is 
used to construct a Pade approximation in the form of a rational function in the parameter 
e. In the third step, the various powers of e which appear in the Pade approximation are 
replaced by new (unknown) parameters {6j}. These new parameters are determined by 
requiring that the residual formed by substituting the new approximation into the governing 
differential equation is orthogonal to each of the perturbation coordinate functions used in 
step one. The technique is applied to model problems involving ordinary or partial differential 
equations. In general, the technique appears to provide good approximations to the solution 
even when the perturbation and Pade approximations fail to do so. The method is discussed 
and topics for future investigations are indicated. 

1 This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-19480 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 


l 




1. Introduction. We have introduced and applied a hybrid perturbation-Galerkin 
technique to a variety of problems involving ordinary differential equations (Andersen 
and Geer, 1991; Geer and Andersen, 1989b, 1990, 1991b), partial differential equations 
(Geer and Andersen, 1991a; Singler and Geer, 1993), and integral equations (Geer and 
Andersen, 1989a), which contain a parameter e. For many different classes of problems, 
this technique appears to significantly improve the usefulness of perturbation expansions 
over a wide range of parameter values. 

In this paper, we show how some of these hybrid ideas can be combined with 
Pade (or rational function) approximations to provide approximate solutions of even 
greater accuracy. The hybrid technique we shall describe successively uses the regular 
perturbation expansion method (see, e.g., Nayfeh, 1973), the Pade expansion method 
(see, e.g. Baker, 1975), and then a Galerkin approximation (see, e.g., Galerkin, 1915). 
In the first step of the method, the regular perturbation method is used to construct 
an approximate solution in the form of a finite power series in a small parameter e 
associated with the problem. In the second step of the method, the series approximation 
obtained in step one is used to construct a Pade approximation in the form of a rational 
function in the parameter e. In the third step, the various powers of e which appear in 
the Pade approximation are replaced by new (unknown) parameters {£_,}. These new 
parameters are determined by requiring that the residual formed by substituting the 
new approximation into the original governing differential equation is orthogonal to 
each of the perturbation coordinate functions used in step one. This hybrid method 
has the potential of overcoming some of the drawbacks of the perturbation method, the 
Pade method, and the Galerkin method when they are applied by themselves, while 
combining some of the good features of each. 

The Pade method has been popular among physicists for improving the conver- 
gence of series or for obtaining information from divergent series (see, e.g., Baker, 1975; 
Gaunt and Guttmann, 1974; and Hunter, 1973). They are easily constructed from the 
first few terms of a power series representation of a function, requiring only the solution 
of a system of linear algebraic equations. From a mathematical point of view, the power 
series representation of a function fails to converge because of the presence of one or 
more (real or complex) singularities of the function, which the series representation is 
unable to adequately approximate. The more general Pade, or rational function, ap- 
proximation of the function can better simulate the presence, or at least the effect, of 
these singularities. This is accomplished through the presence of zeros in the denomi- 
nator of the approximation (the simulation of poles) or the combination of zeros of the 
numerator and denominator at nearly the same locations (which can simulate certain 
branch cuts) (see, e.g., Andersen and Geer, 1982; Dadfar, et al, 1984; or Van Dyke, 
1974). However, the region of convergence of a sequence of Pade approximants is, in 
general, not closely related to the region of convergence of the power series on which 
the approximations are based. In fact, the construction of Pade approximants in the 
usual way may lead to approximations which fail to converge at points where the Taylor 
series of the function does converge. Often this is due to the appearance of spurious 
poles or zeros in the approximants which are not present in the original function. (See, 
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e.g., Bender and Orszag, 1978, for some simple, but very interesting, examples which 
illustrate this phenomena.) As we shall demonstrate below, the hybrid method we will 
describe appears to overcome this drawback of the Pade method by “adjusting” the lo- 
cation of the poles and zeros of the approximants to simulate better the true analytical 
structure of the function being modeled. 

Some general observations about the hybrid technique we are proposing are the 
following. First, in many perturbation problems, much effort has to be expended to 
compute each additional term in a perturbation expansion. Through the use of the 
proposed hybrid method, the “information” contained in the known perturbation terms 
can be exploited more fully. In particular, more of the analytical structure of the 
solution can be “uncovered” using a rational function approximation. Secondly, another 
way of viewing the technique is to recognize that in many perturbation expansions 
the functional form of the higher order terms can be well approximated by a linear 
combination of the lower order terms. Thus, much of the effect of the higher order 
terms may be included by applying the technique to the lower' order terms. Finally, 
preliminary unpublished investigations indicate that, while the use of a Taylor series 
expansion is frequently limited by a finite radius of convergence, the proposed hybrid 
method can sometimes yield good results even well outside the radius of convergence. 

In Sec 2 we shall describe our method in more detail and then apply it to two 
model boundary value problems involving ordinary differential equations (Sec 3) and 
to two model problems involving partial differential equations (Sec 4). Finally, we shall 
discuss observations about the method in Sec 5, and also indicate topics for further 
investigation. 

2. Description of the Method. We consider the problem of finding approxima- 
tions to the solution u(x,e) to the problem 

(1) L(u , e) = 0, x G D, with B(u, e) = 0, x 6 dD. 

Here L represents a (linear or nonlinear) differential operator, x is a scalar or a vector, 
D is the domain in which the differential equation is to hold, dD denotes the boundary 
of D, and e is a small parameter. For simplicity we shall assume that the operator B , 
which defines the boundary conditions of the problem, is linear and homogeneous in u. 

The hybrid Pade-Galerkin technique we wish to present can be conveniently de- 
scribed as the following three-step process. 

Step one: For small values of the parameter e, we assume that u can be expanded 
in a (convergent or asymptotic) perturbation series of the form 

(2) u{x, e) = E u i (x)e J + 0(e n+1 ). 

i=o 

Here each of the perturbation coordinate functions Uj is determined in a straightforward 
manner using the regular perturbation method. 

Step two : We now use the functions {u ; } determined in step one to construct the 
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Pade approximant P[M , N] defined by 


( 3 ) 


P[M,N J(x,e) 




w 0 (x) = 1. 


Here the functions {u?(a:)} and {u>j(x)} are determined by the condition that P[M, iV](m, c) 
u(x, e) -f 0(e M+N+1 ), as e — *■ 0. This condition leads to the requirements 


( 4 ) 


Uj(x) — j — 1) •••> 

k = 0 


min(p,JV) 

(5) ^2 u P-k(^)wk{x) = -itp(x), p = M + 1, M 4- 2,...,M + N. 

k = 1 


Step three: Once P[M, iV] has been determined, we define the hybrid Pade-Galerkin 
approximation H [M, iV] by 


( 6 ) 


H[M,N](x,e) 


EjLo 

1 + Ef=, wj(x)S MH ' 


Here H[M,N] has the same functional form as P[M,N], except that each power of 
e in the definition of P[M , N] (except the coefficient of wq) has been replaced by a 
new (unknown) parameter, or amplitude, 6j = <5j(e). The M + N + 1 amplitudes {(5j} 
are determined by requiring that the residual formed by substituting H [M, N] into 
the governing differential equation is orthogonal to each of the perturbation coordinate 
functions used in step one. Thus, the new amplitudes {<5j} are determined from the 
Galerkin conditions 

(7) / L ( H[M , N]{x , e), e) u k dx = 0, k = 0, 1, M + N. 

Jd 

Conditions (7) are a system of A/ + N + 1 equations for the M + N + 1 amplitudes {<5j}. 
In general, these equations are nonlinear and must be solved numerically. However, 
we note that they can be solved efficiently using Newton’s method, starting with small 
values of e, where we expect 6j = e J , for 0 < j < M, and = e J , for 1 < j < N, 
and then incrementally proceeding to larger values of e. 

Before applying this technique, we make a few observations. In step one of the 
method, i.e., the perturbation step, the solution to a particular problem involving a 
small parameter is developed in terms of a series of unknown functions with preassigned 
coefficients, i.e. gauge functions. The unknown functions are usually determined by 
solving a recursive set of differential equations which are, in general, simpler than the 
original governing differential equation. Using these perturbation coordinate functions, 
a new approximation in the form of a rational function of e can be constructed in a 
straightforward manner. (We note that the Pade approximant P[M, 0] as defined above 
is the same as the perturbation approximation of order M, while P[M, iV] with N > 0 is 
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a truly new approximation.) By contrast, using the Galerkin technique in the standard 
manner, one seeks an approximate solution to the problem in the form of a series of 
specified (known) coordinate functions with unknown coefficients. The coefficients are 
determined by requiring that the residual formed by substituting the trial solution into 
the governing differential equation is orthogonal to each of the coordinate functions. The 
technique we are proposing simply replaces the power series form of the approximation 
in the Galerkin step by a more general rational function approximation to the solution. 

We also observe that the perturbation, Pade, and Galerkin methods have each been 
useful when applied by themselves in providing approximate solutions to a wide variety 
of nonlinear (and otherwise difficult) problems. However, each of these techniques has 
certain drawbacks. The perturbation method has at least two major drawbacks. First, 
as the number of terms in the perturbation expansion increases, the mathematical 
complexity of the equations which determine the unknown functions increases rapidly. 
Thus, in most practical applications, the perturbation series is limited to only a few 
terms. A second drawback to the perturbation method is the requirement of restricting 
the perturbation parameter to small values in order to obtain solutions of acceptable 
accuracy. (These drawbacks of the perturbation method have been recognized and 
several modifications or extensions have been proposed, as, e.g., in Andersen and Geer, 
1982, and Van Dyke, 1974.) One of the main drawbacks of the Pade method is that it 
often has to be applied ” blindly”, with little or no assurance that the approximations 
obtained will converge to the solution as the number of terms is increased. Also, a 
rational function obviously has (real or complex) poles which may or may not improve 
the quality of the approximations obtained. The main shortcoming of the Galerkin 
method is the difficulty, from a practical point of view, of selecting a small number of 
good coordinate functions. 

By basing our technique on a formal perturbation expansion of the solution, we 
should note that our particular choice of coordinate functions overcomes the main draw- 
back of the Galerkin method (at least for the case when N = 0.) By the way they are 
constructed, the perturbation coordinate functions are (under certain assumptions) el- 
ements of a set of functions which span the space of solutions in a neighborhood of 
their point of generation. Thus, they should fully characterize the solution u in that 
neighborhood. Also, in many applications, the functions {ttj} are determined by solving 
a set of linear equations, even though the original operator L may be nonlinear. The 
first property is necessary for the convergence of the Galerkin method, while the second 
property enhances the effectiveness of the proposed hybrid method for solving nonlinear 
problems. In addition, for TV > 0, we shall demonstrate below that the Galerkin method 
appears to determine the parameters of the Pade approximants in such a manner that 
the poles of H [Af, N] tend to improve and not degrade the quality of the approximations 
obtained. 

In the following sections we shall apply our method to several model boundary value 
problems involving ordinary or partial differential equations. These examples will serve 
to illustrate the application of the method, and will also provide insights into some of 
the mathematical properties of the method. A more thorough and detailed analysis of 
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some of the mathematical properties of the technique, as well as the application of the 
method to certain singular perturbation problems, will be presented elsewhere. 

3. Applications to Ordinary Differential Equations. In this section we shall 
apply our method to two model boundary value problems involving second order ordi- 
nary differential equations. 

Example 1 : We consider first the two-point boundary value problem 

(8) L(u , e) = u"(x) — eu'(x) — e = 0, 0 < x < 1, 


with B(u, e) = u = 0, for x = 0 and x = 1. 

To help evaluate the accuracy of our proposed hybrid method, we note that this simple 
problem has the exact solution 

1 - e £X 

(9 ) u =~r^- x - 

The Taylor series expansion of u(x,e) about e= 0 converges for |e| < 27T, since the 
singularities of u which lie closest to the origin in the complex e— plane are located at 
e = ± 2ni. We also note that u develops a boundary layer at x = 1 as e becomes large. 

Step one: Using the regular perturbation method we find that each of the perturba- 
tion coefficient functions Uj in Eq (2) is a polynomial of degree j + 1 which vanishes at 
x = 0 and x = 1. In particular, we find u 0 = 0, ui = x(x— 1)/2, i <2 = x(x — l)(2x— 1)/12, 
u 3 = x 2 (l — x) 2 /24, and u 4 = x(x — l)(2x — l)(3x 2 — 3x — l)/720. 

Step two: We now use the perturbation coordinate functions {uj(x)} determined 
in step one to calculate Pade approximations P[M, N] for various values of M and N. 
In particular for M + N = 2, using Eqs (4) and (5), we find 

(10) P[ 2,0] = eui(x) + e 2 u 2 (x), 


which is the second order perturbation solution, and 


( 11 ) 


P[l,l] = 


ex(x — l)/2 
1 + e(l — 2x)/6’ 


since Vo = 0, v\ = x(x — 1 )/6, and w\ = (1 — 2x)/6. We note that P[l, 1] has a pole 
at x = 1/2 + 3/e, which lies inside the interval [0, 1] for e > 6. Thus, we do not expect 
P[l, 1] to be a good approximation to u(x, e) for e> 6. In Figs 1 and 2, we have plotted 
P[2, 0] and P[l,l], along with the exact solution, for e = 5 and e = 10, respectively. 
The influence of the pole of P[l, 1] is evident in the Figure 2. 

Step three: We now use the form of the Pade approximations determined in step 
two to define the hybrid approximations 


( 12 ) 


H[ 2, 0] = <5iui(x) + 6 2 u 2 {x), 
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and 


(13) 


tt r-i ii _ <^l x ( x \)J2l 

l ’ J 1 + <5 2 (1 — 2x)/6’ 


where the amplitudes 8 \ = ^i(e) and 82 = 82 (c) are determined from the Galerkin 
conditions (7), which for this example become 

(14) f 1 ( H"[M,N } - eH'[M,N } - e)u k dx = 0, for k = 1,2. 

Jo 

When M = 2 and N = 0 Eqs (14) are the linear equations 

(1/12)$! + (c/720)8 2 = e/12, 

(e/720)$i — (l/720)$2 = 0, 


from which we find 
(15) 


C _ C C = L 

1 1 + e 2 /60’ 2 1 + e 2 /60 


When M = 1 and N = 1, from Eqs (14) we find that 8 1 and <$ 2 satisfy the nonlinear 
equations 


(16) 


2e$| 

9(3e + $ 2 )(36 - 81) log (|±|) - 324e<5 2 - 108$| + 6e^’ 


(17) 


/6 + $ 2 \ 4$ 2 (324e + 216$ 2 - 9e8j - 4 $|) 

V6 -$ 2 / (36 — $2)(108e + 725 2 — e8%) 


Equation (16) expresses 8 1 in terms of 8 2 , while Eq (17) is an equation for 8 2 as a 
function of e, which must be solved numerically. In Figs 1 and 2 we have also plotted 
H[2, 0] and J/[l,l] for e = 5 and e = 10, respectively. As these figures illustrate, the 
hybrid approximations are more accurate than the Pade approximations on which they 
are based. 

Although we shall discuss, more fully, insights about our method provided by this 
example (as well as our other examples) in Sec 5, we make a few preliminary observations 
at this point. First, from Eqs (15) and Eqs (16)-(17) it follows that 


81 = e + 0(e 3 ) and 8 2 = e 2 + 0(e 3 ), as e — ► 0, for M = 2 and N = 0, 


(18) $1 = e + 0(e 2 ) and 8 2 = c + 0(e 2 ), as e — > 0, for M = 1 and N = 1. 

Thus if[2,0] reduces to P[2,0] and #[1,1] reduces to P[l, 1], as c — > 0. In addition, as 
e — > 00, we find from Eqs (16) and (17) that, for the case M = 1 and N = 1, 

$i=4 + 0( 7(e)) and 8 2 — 6 - 7(e) + o(7(e)), 
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as £ -> 00. 



Here 7(e) -> 0 + as e -> +00 and is defined by 7log(7) = -12/e. From these results, as 
well as the numerical solution of Eqs (16) and (17) for “intermediate” values of e, we 
find that the pole of H[ 1, 1], which is located at x = 1/2 + 3/^, lies outside the interval 
[0,1] for all positive values of e. (The location of the pole approaches 1 from above as e 
becomes large, with the pole lying, for example, at about 1.008 when e = 50.) 

It is straightforward to use the procedure and formulas outlined above to construct 
Pade and hybrid approximations to u(x , e) for other (larger) values of M and N. For 
example, we find 


P[ 4, 0](x, e) = Oii(x) + e 2 u 2 (x) + e 3 u 3 (x) + e 4 u 4 (x), 


H[ 4, 0](x, e) = ^iUi(x) + 6 2 u 2 {x) + S 3 u 3 (x) + 6 4 u 4 (x), 


where 


15120c + 420e 3 
15120 + 420e 2 + e 4 ’ 


8 2 = e<5i, S 3 


15120c 3 

15120 + 420e 2 + e 4 ’ 


8 4 — eS 3 ] 


and 


P[ 2,2] 


evijx) + e 2 u 2 (x) . . _ ^iui(g) + ^2^2(3^) 

1 + etyi(x) + ’ 1 + 8 3 wi(x) + 6 4 w 2 (x)’ 


with 


v\ = x(x — l)/2, v 2 


x(x — l)(2x — l)(x + l)(^r — 2) 
60(x 2 — x + 1) ’ 


Wl 


(2x - l)(2x 2 - 2x + 1) 
10(x 2 — x + 1) 


and 


w 2 


3x 4 — 6x 3 + 4x 2 — x + 1 
x 2 — x + 1 


In Fig 3 we have plotted the Pade approximation P[2, 2], the hybrid approximations 
H[ 4,0] and H[ 2,2], as well as the exact solution, for e = 50. Here H[ 2,2] is clearly the 
“best” of the approximations shown and does a good job of simulating the boundary 
layer behavior of the exact solution, even at this moderately large value of e where the 
boundary layer is well formed. 

Example 2: We consider now the two-point boundary value problem 


(19) Z(u, e) = u"(x) + e sin(2x)i/(x) + 2c cos(2x)u — 2e cos(2x) = 0, 0 < x < 7r, 


with B(u , e) = u = 0, for x = 0 and x = 7 r. 
The exact solution to this problem is 
(20) u = 1 - e - £8in2(x) . 
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The Taylor series expansion of u(x, e) about e= 0 converges for all |e| < oo, since u has 
no singularities in the finite part of the complex e— plane. For this example, u develops 
boundary layers at both x = 0 and x = 7r as e becomes large. 

Step one: Using the regular perturbation method we find the Uj in Eq (2) are given 
by u 0 (x) = 0 and uj(x) = (— 1) J+1 sin 2j (a;)/j!, for j > 1. 

Step two: Using the perturbation coordinate functions (uj(x)} determined in step 
one we can calculate Pade approximations P[M, iV] for various values of M and N. In 
particular for M + N = 2, using Eqs (4) and (5), we find 

(21) P[ 2, 0](x, e) = esin 2 (x) — e 2 sin 4 (x)/2, 

which is the second order perturbation solution, and 


( 22 ) 


P[l,l](x,e) 


esin 2 (x) 

1 + (e/2)sin 2 (x)’ 


We note that, for t > 0, the poles of P[l,l] are all complex and, hence, lie outside 
the interval D = [0, 1]. In particular, the poles which lie closest to D are at x = 
±i sinh -1 (-y/2/e) and x = 7r ± i sinh _1 (y^/e), which approach x = 0 and x = ir, 
respectively, as e -» +oo. In Fig 4, we have plotted P[2,0] and P[l, 1], along with the 
exact solution, for e = 10. The influence of the (complex) poles of P[l, 1] in helping to 
simulate the steep slope of the solution near the ends of the interval is evident in Fig 4. 

Step three: We now use the form of the Pade approximations determined in step 
two to define the hybrid approximations 


(23) 
and 

(24) 


H[ 2,0] = 6isin 2 (x) — £2 sin 4 (x)/2, 


mi i] (* c ) 

1 ’ J_ 1 + (^2/2) sin 2 (x) ’ 


where the amplitudes 61 and 62 are determined from the Galerkin conditions (7), which 
for this example become 

/ 1 {H[M, N}" + esin(2 x)H[M, N]' + 2ecos(2 x)H[M, N] 

Jo 


(25) — 2ecos(2x)} Ukdx = 0, for k = 1,2. 

When M = 2 and N = 0 Eqs (25) are linear and their solution is 

(26) = 16 i( ? + c ) - 6 = . 

v ' 32 + 16e + 3e 2 2 32 + 16e + 3e 2 

When M = 1 and N = 1, Eqs (25) are nonlinear and must be solved numerically. In 
Fig 4 we have also plotted /7[2,0] and H[ 1, 1] for t — 10. As the figure illustrates, the 
hybrid approximations are again more accurate than the Pade approximations on which 
they are based. Using Eqs (25), we can again show that the relations (18) hold for this 
example as well. Also, it is straightforward to use the procedure and formulas outlined 
above to construct Pade and hybrid approximations to u(x, e) for other (larger) values 
of M and N. 
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4. Applications to Partial Differential Equations. In this section, we shall 
apply our method to two model elliptic boundary value problems in two space dimen- 
sions. (Although we shall limit our illustrations to these two dimensional examples, the 
method itself is formally applicable to problems involving any number of dimensions.) 

Example 3: We consider the elliptic boundary value problem 

L(u , e) = V 2 u + ecos(x) sin^it* + esin(x) cos (y)u y 


(27) 


— 2e sin(x) sin (y)u + 2e sin(x) sin(y) = 0, 


for (x, y) e D = {(x, y ) : 0 < x < tt, 0 < y < 7r}, with B(u, e) = u = 0 for (x, y) € dD. 
Here V 2 denotes the usual (two dimensional) Laplacian operator. This problem has the 
exact solution 


(28) 


U = 1 — e ~ « sin(ac) slti(v) ^ 


and hence its Taylor series expansion about e = 0 converges for all |e| < oo. We also 
note that, as e — ► +oo, u develops a boundary layer around the entire boundary of D , 
and approaches the value of 1 at each point in the interior of D. Figure 5 depicts a 
surface plot of the exact solution for e = 10. 

Step one : Using the regular perturbation method we find 


(29) 


u o 


= 0, and Uj = - — r . — sin J (x) sin-'(y), for j > 1. 


Step two: Using the perturbation coordinate functions in Eqs (29) we find 


(30) 


P[ 2, 0] = e sin (x) sin (y) — — sin 2 (x) sin 2 (y), 


(31) 


PM 


esin(x)sin(y) 

1 + (e/2) sin (x)sin(y)' 


We note that P[l,l] has real poles for e > 2 where sin (x) sin (y) = —2/e, and hence 
these poles lie outside of D for all values of e> 0. In Fig 6 we have plotted a cross- 
section of the approximate solutions P[2, 0] and P[l, 1], along with the exact solution, 
at y = 7t/2 for 0 < x < tt for e = 10. 

Step three : Using the form of the Pade approximants above, we define 


(32) 


H [2, 0] = di sin (x) sin (y) — 


. 

— sin 

2 


2 (x)sia 2 (j/), 


H{ 1,1] 


di sin (x) sin (y) 

1 + (8 2 /2) sin (x) sin (y) ’ 
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(33) 



where 8\ and S 2 are determined from the Galerkin conditions (7), which for this example 
become 


(34) f f L(H[M, N],e)uk(x,y)dxdy = 0, fork = 1,2. 

Jo Jo 

For N = 2 and TV = 0, these equations are linear and yield the solutions 

e 32e(6075;r 4 + 7632;r 2 e - 409600) c 800e 2 (243;r 4 - 16384) 
S ' A ’ *' = A ’ 


(35) 


A = 194400tt 4 - 13107200 + 244224;r 2 e + (2097152- 18225rr 4 )e 2 . 


From Eqs (34) and (35) it follows that 8\ and S 2 satisfy the conditions (18) as e — » 0. 
In Fig 6 we have also plotted the approximate solutions H[2, 0] and #[1,1] at y = 
7t/ 2 for 0 < x < 7r. The improvement of the hybrid solutions over the corresponding 
perturbation solutions is evident from the figure. 

Example 4: We consider the problem 


(36) 


L(u , e) = V 2 u + e sin(x) cos (y)u x + 2 e sin(x) sin (y) = 0, 


for ( x,y ) € D = {( x,y ) : 0<x<7r, 0 < 2 / < x), with B(u,e) = u = 0 for 
( x ,y) G dD. Although there appears to be no closed form solution for this problem, 
it is straightforward to show that its solution u(x , y ) exhibits a number of interesting 
properties. The solution is positive over the entire domain D for positive values of e 
and is invariant under the 180° rotation x — > w — x, y — ► n — y. Since Eq (36) is also 
invariant under the transformation e — ► — e, and u(x,y) — » —u(x,tt — y ), we may focus 
our attention on solutions for positive e. Figure 7 depicts a surface plot of the “exact” 
solution, obtained by purely numerical means, for e = 12. The solution has a regular 
perturbation expansion about t = 0, which converges for values of t with |e| < e 0 — 6.7. 
(See Geer and Andersen, 1991a, where this equation was discussed in some detail.) In 
addition, ase-+ + 00 , the solution develops a number of boundary layers over portions 
of dD , but not over the entire boundary. 

Step one: Using the regular perturbation method, we find 

tt 0 = 0, ux = sin(x)sin(y), u 2 = (l/32)sin(2x) sin(2y), 

U 3 = Y^g |^sin(3x) sin(3y) + | sin(3z) sin(y) - ^ sin(ar) sin(3y) - sin(a:) sin(y)| , 


u 4 = 


1 


49152 


sin(4x) sin(4y) + ^ sin(4x) sin(2y) - 


1 


19200 


sin(2x)sin(4y) 


1920 


sin(2x) sin(2y). 


(37) 
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Step two: Using Eqs (37) we find 


P[ 2 , 0 ] = esin(x) sin(y) + — sin( 2 x) sin( 2 y), 


( 38 ) 


P[hl] 


esin(x) sin (y) 

1 — (e/ 8 ) cos(x) cos(y) ’ 


We note that the poles of P[l,l] lie inside D for e > 8 . In Fig 8 we have plotted 
approximations to the diagonal cross-section y = x of the solution at e = 12 using 
P[2 , 0 ], P[ 1, 1 ], and the numerical solution. 

Step three: Using the expressions in (38) we define 


(39) 


£ 

H[ 2 , 0 ] = (5i sin(x) sin (y) + sin( 2 x) sin( 2 y), 


sin(x) sin(y) 

’ 1 — (^ 2 / 8 ) cos(x) cos(y)’ 


For M = 2 and N = 0 , the Galerkin conditions (7) are linear and give 

A 


(40) 


6x = 


r. *2 = 


1 + e 2 /128’ 1 + e 2 /128’ 


For M = 1 and JV = 1 , these conditions are nonlinear and must be solved numerically. 
Using these numerically computed solutions, it appears that the poles of H[ 1 , 1 ] always 
lie outside of D. In Fig 8 we have also plotted approximations to the diagonal cross- 
section y = x of the solution at e = 12 using #[2,0] and H[ 1, 1]. 

The quality of the hybrid approximations improves as the order of the Pade ap- 
proximants is increased. In particular, we find 


and 


where 


p [ 4 > 0] = S y)> H \ 4 > °] = Yj y)> 

j = 1 j = 1 


pro 91 - et?1 + ^ rrro 9 1 _ Vl + 
L,J 1 + e Wl + e>w 2 ' l,J + 


wi = sin(x) sin(y), 


(65 + cos(2x) + 193cos(2y) — 25cos(2x) cos(2t/)) sin(2x) sin(2y) 
160(125 — 11 cos(2x) + 61 cos(2y) + 5 cos(2x) cos(2y) 
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W 1 = 


1201 cos(x) cos(y) — 31 cos(3x) cos(y) + 137 cos(x) cos(3y) + 25 cos(3x) cos(3y) 
80(125 — 11 cos(2x) + 61 cos(2y) + 5 cos(2x) cos(2y) 


W 2 = {32301 — 7590cos(2x) + 237 cos(4x) + 25962cos(2y) — 4396cos(2x) cos(2y) 
+106cos(4x) cos(2y) + 2109cos(4y) — 902cos(2x) cos(4y) + 125cos(4x) cos(4y)} / 

(23040(125 — 11 cos(2x) + 61 cos(2y) + 5 cos(2x) cos(2y)). 

The amplitudes {<5j} are determined as indicated above. In particular, for M = 4 and 
N = 0 we find 

56524800c + 1216128c 3 

1 “ 56524800 + 1216128c 2 + 2141c 4 ’ 2 ~ 6 u 

56524800c 3 _ 

3 “ 56524800 + 1216128c 2 + 2141c 4 ’ 4 ~ c 3 ’ 

while the amplitudes for M = N = 2 must be determined numerically. 

In Fig 9 we have plotted approximations to the diagonal cross-section y = x of the 
solution at e = 12 using P[ 4, 0], H[ 4, 0], P[ 2, 2], and H[ 2, 2]. From the figure we see that 
each hybrid approximation is an improvement over the Pade approximation on which 
it is based, while H [2, 2] appears to be the “best” of all of the hybrid approximations. 

5. Discussion and Observations. We now make some observations concerning 
our method which have been provided by the examples we have presented. 

In Ex 1 we considered a simple two point boundary value problem. The hybrid 
solutions presented provide approximations which are consistently more accurate than 
the perturbation or Pade approximations on which they are based, and they appear to 
be converging to the exact solution as the number of terms is increased. This is true 
for small values of e, for which the perturbation and Pade solutions give reasonable 
approximations to the solution, as well as for larger values of e, where the perturbation 
and Pade solutions fail to give reasonable approximations. Part of the reason for this 
behavior is related to the convergence properties of the original perturbation series. The 
Taylor series expansion (and, hence, the regular perturbation expansion) of the solution 
converges only for values of e with |e| < 2n. Thus, the perturbation solution, by itself, is 
useless for computing approximations to u(x, e) for values of e with |e| > 2n. However, we 
have shown (Andersen and Geer, 1991) that the sequence of approximations {H[M, 0]} 
converges to the exact solution as M — ► oo for each 0 < x < 1 and for each 0 < e < oo. 
In addition, for N > l, each Pade approximation P[M,N] fails to give a reasonable 
approximation to the solution when one or more of its poles moves inside the domain 
(interval) D = (0, 1). The presence of a pole in a Pade approximation is undoubtedly 
related to the formation of the boundary layer at x = 1 as e increases. When this pole 
lies outside of D , but “close to” x = 1, it helps to simulate the steep slope of the solution 
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forming within the boundary layer. When a pole moves inside D, however, it creates a 
singularity in the approximation which is not present in the exact solution. It appears 
that the poles of the hybrid approximations H [M, iV] always lie outside of D and hence 
H [M, N] has the potential of providing at least a reasonable approximation to the exact 
solution, even when the corresponding Pade approximation cannot do so. In Fig 10 we 
have plotted the location of the pole of P[ 1 , 1] and of H[ 1 , 1] for 0 < e < 30. As the 
figure illustrates, the pole of H[ 1, 1] asymptotically approaches 1 from outside of D, as 
e — *• oo. 

In Ex 2, the perturbation solutions converge to the exact solution for all finite values 
of e, as the number of terms in the approximation is increased. However, as e increases, 
boundary layers form at each end of the domain (interval) D, and, consequently, more 
terms are required in the perturbation expansion to provide approximations of “accept- 
able” accuracy. In this case, the poles of the Pade and hybrid approximations lie in 
the complex plane and hence outside of D for all real values of e. Some of these poles 
do approach the ends of the interval D = [0, 1], however, as e -> oo, and hence help 
to simulate the actual behavior of the exact solutions in these regions. In Fig 11 we 
have indicated the location of these poles of P[l, 1] and H[ 1, 1] as a function of e. From 
the figure it is clear that the poles of H[ 1, 1] lie closer to the boundary of D than do 
the corresponding poles of P[l, 1], This helps to explain the increased accuracy of the 
hybrid approximations. 

In our third example, we examined a two dimensional, elliptic boundary value 
problem involving a simple partial differential equation. As in Ex 2, the perturbation 
solution converges to the exact solution for all finite values of e, as the number of terms 
in the approximation is increased. However, as e increases, boundary layers now form 
over the entire boundary of the domain D. Consequently, more terms in the perturbation 
expansion are again required to provide approximations of “acceptable” accuracy as t 
increases. In this case, all of the real poles of the Pade and hybrid approximations 
P[l, 1] and H[ 1, 1] lie outside of D for all real values of e. In Fig 12 we have plotted the 
location of the real poles of these approximations for selected values of e > 2. From the 
figure we see that, as in Ex 2, the poles of H[ 1 , 1] lie closer to the boundary of D than 
do the corresponding poles of P[l, 1]. 

In Ex 4 we considered another two dimensional, elliptic boundary value problem, 
whose exact solution develops a more complicated boundary layer structure as e in- 
creases. In this case the perturbation solution has a finite radius of convergence and 
hence the classical perturbation approximations are useless for values of e with magni- 
tude greater than the radius of convergence. However, all of the hybrid solutions appear 
to give reasonable approximations even for values of e moderately above the radius of 
convergence (see Figs 8 and 9). In this case, the poles of P[l,l] move inside D for 
values of e > 8. The location of these poles are indicated in Fig 13. The corresponding 
poles of H[ 1, 1] are purely imaginary (since |6 2 | < 8) and hence they lie outside of D. 

The observations we have indicated above, as well as experience we have had with 
the method for other examples, motivate us to make the following conjecture. It appears 
that, if the original perturbation expansion has a finite radius of convergence, then the 
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poles of the Pade approximations constructed from this expansion will move inside 
the domain D as e increases. If the original perturbation expansion has an infinite 
radius of convergence, then the poles of the Pade approximations will remain outside 
of D. In the first case, the hybrid approximations we have described appear to keep 
the corresponding poles outside of D for all values of e, while in the second case the 
hybrid solutions appear to move the poles closer to the boundary of D. In either case, 
the hybrid technique appears to “adjust” the location of the poles in such a manner 
as to provide approximate solutions of greater accuracy than the corresponding Pade 
approximation. Obviously, this conjecture needs to be studied in more detail and is the 
subject of some of our current investigations. 

6. Conclusions. In general, the basic idea of using the coefficient functions in the 
original perturbation expansion in a small parameter e to construct one or more Pade 
approximations seems to be justified, and, in fact, is recommended when the solution 
tends to exhibit some kind of boundary layer behavior as e increases. The increased 
ability of the Pade solutions to simulate boundary layer behavior appears to provide 
more accurate approximations, even for “smaller” values of e, for which the boundary 
layer is not yet well formed (see, e.g., Figs 1 and 2). For “any” value of e, however, 
the hybrid solutions are consistently more accurate than the Pade solutions on which 
they are based and it is recommended that the additional step of computing the new 
amplitudes {^} be undertaken, especially when only a few perturbation coordinated 
functions are known. 
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0 x 1 


FIG 1 . The Pade approximations P[2, 0] and P[l,l], the hybrid approximations 
J/[2, 0] and H[ 1, 1], and the exact solution (circles) for Ex 1 with e = 5. 



I 1 1 1 i i i» i i 1 1 

0x1 


FIG 2. The Pade approximations P[2,0] and P[l, 1], the hybrid approximations 
H[ 2,0] and f/jl, 1], and the exact solution (circles) for Ex 1 with e = 10. For this value 
of e, the pole of P[l, 1] lies at x = 0.8. 
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FIG 3. The Pade approximation P[ 2,2], the hybrid approximations #[4,0] and 
H[ 2,2], and the exact solution (circles) for Ex 1 with e = 50. 



FIG 4. The Pade approximations P[2, 0] and P[l,l], the hybrid approximations 
H[2, 0] and H[ 1, 1], and the exact solution (circles) for Ex 2 with t — 10. 
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FIG 5. A surface plot of the exact solution (Eq (28)) for Ex 3 with t = 10. Note 
the rather steep boundary layer around the entire boundary of D. 



FIG 6. The Pade approximations P[2, 0] and P[l,l], the hybrid approximations 
H[2, 0] and If [1,1], and the exact solution (circles) for Ex 3 along the cross section 
y = 7t/2 with e = 10. Note the very “poor” quality of the perturbation solution P[2,0], 
on which all of the other (much better!) approximations are ultimately based. 
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FIG 7. A surface plot of the “exact” solution for Ex 3 obtained by purely numerical 
means for e = 12. Note the boundary layers beginning to form around portions of the 
boundary of D, especially near (x,y) = (0,0) and (x,y) = (7r,7r), as well as an internal 
“rise” forming near the diagonal x = y. 



FIG 8. The Pade approximations P[ 2,0] and P[l,l], the hybrid approximations 
ff[2,0] and J/[l,l], and the exact solution (circles) for Ex 4 along the diagonal cross 
section x — y with t = 12. 
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FIG 9. The Pade approximations P[4, 0] and P[ 2,2], the hybrid approximations 
#[4, 0] and H[ 2,2], and the exact solution (circles) for Ex 4 along the diagonal cross 
section x = y with e = 12. 



0 e 30 

FIG 10. A plot of the location of the poles (vertical axis) of P[l, 1] and H[ 1, 1] 
for Ex 1 as e (horizontal axis) varies between 0 and 30. The “critical” value of x = 1 is 
indicated, as well as the asymptotic value (short dashed line) of the pole of P[l, 1] as 




FIG 11. A plot of the magnitude (vertical axis) of the poles of P[ 1, 1] and H[ 1, 1] 
for Ex 2 which lie closest to x = 0 as e (horizontal axis) varies between 0 and 20. These 
poles are purely imaginary and lie at x = sinh -1 (y2/e) and x = sinh (y 2 /^ 2 ) 
for P[l, 1] and H[ 1, 1], respectively. 



FIG 12. A plot of the locations of some of the poles of P[l, 1] (dashed lines) an <l 
H[l, 1] (solid lines) in the s,y-plane which lie closest to D for Ex 3 wit e = . , , 
and 10. The boundary of D is also indicated. Note that the poles of H[ , ] 1 
to the boundary of D than do the corresponding poles of P[ 1, 1J- 


21 





FIG 13. A plot of the locations of some of the poles of P[l, 1] (dashed lines) in 
the x,y — plane which lie closest to D for Ex 4 with e = 10, 15, and 20. The boundary 
of D is also indicated. The influence of the pole of Pfl, 1] at e = 12 is also evident in 
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